Learning plan

1. The topics in this deck mirror the workflow of a project

flowchart LR
  D[Data exported\nfrom REDCAP]-->C[Data cleaning]
  C-->S[Statistical analysis]
  S-->P[Presentation\nof results]

2. The game plan is to first do a quick pass through from data cleaning to analysis:
   1. Take a toy dataset
   2. “Clean” the data (minimally)
   3. Run a regression model
3. Once we have gotten comfortable with R, we will explore each topic in greater detail
   1. Data cleaning
   2. Programming fundamentals
   3. Analysis: generalised linear models
   4. Presentation of results

What is an R package?

1. It is a collection of code/data written by someone and packaged for distribution
2. Hosted on the Comprehensive R Archive Network and/or Github for public download
3. From the CRAN webpage, we are able to find a reference manual documenting the data/functions of the package
4. From the Github repository, we are able to see the actual codes implemented by the package’s author(s)
5. To install a package we would run install.packages("palmerpenguins") in the console
6. To load a package, we would run library(palmerpenguins)
7. To uninstall a package we would run remove.packages("palmerpenguins") in the console

[DEMO] Intro to RStudio workspace

1. Console: where you run/execute lines of code

2. To assign a variable, we use the <- operator

   x <- 1:10; print(x)

    [1]  1  2  3  4  5  6  7  8  9 10

   1. random-access memory is allocated when a variable is assigned/declared
   2. random-access memory is freed-up when the variable is deleted (or garbage collected if it goes out of scope)

3. Environment: where you’re able to see variables stored in random-access memory

   1. To clear variables from environment \(\rightarrow\) execute rm(list = ls()) in console

[DEMO] Intro to RStudio workspace

1. Script
   1. Writing code and saving the script
   2. Execute single/multiple lines of code \(\rightarrow\) highlight code + Ctrl-enter
   3. Runs from top to bottom
   4. Comments will not be executed
2. Hotkeys
   1. Clear the console \(\rightarrow\) Ctrl-l
   2. Clear a line in console \(\rightarrow\) Ctrl-d

[PRACTISE] Tasks to do

1. For this task, you will need the palmerpenguins and tidyverse packages. Install them if you have not already done so
2. We will make use of the penguins dataset from palmerpenguins package
3. Understand the data
4. Dichotomise body_mass_g (continuous) into a categorical variable with 2 categories (“light”: <= 4750g and “heavy”: > 4750g)
5. Run a linear regression¹ with body_mass_g as the dependent variable and independent variables: species, bill_length_mm, bill_depth_mm, flipper_length_mm and sex
6. Run a logistic regression with the categorical body mass variable you created earlier as the dependent variable, with the same independent variables as before

Data types

Data types (continued)

1. Factor datatype for working with categorical variables

   x <- factor(sample(c("Malay", "Others", "Indian", "Chinese"), 10, replace = T),
               levels = c("Chinese", "Malay", "Indian", "Others"))
   print(x)

    [1] Indian  Malay   Malay   Chinese Chinese Chinese Malay   Malay   Others 
   [10] Others 
   Levels: Chinese Malay Indian Others

   print(levels(x))

   [1] "Chinese" "Malay"   "Indian"  "Others" 

2. A factor variable is more useful than a plain character vector. E.g. the first level will be used by glm as the reference category

Data structures (vectors)

Data structures (lists)

1. Vectors can only hold 1 kind of data type e.g. characters, numeric etc
2. Lists are “containers” that can hold multiple data types

   1. Unnamed lists

      list("a", 1, 10:15)

      [[1]]
      [1] "a"
      
      [[2]]
      [1] 1
      
      [[3]]
      [1] 10 11 12 13 14 15

   2. Named lists

      list(a = "a", b = 1, c = 10:15)

      $a
      [1] "a"
      
      $b
      [1] 1
      
      $c
      [1] 10 11 12 13 14 15

Data structures (dataframe)

Data structures (dataframe, accessing variables)

1. To get the column/variable-names of your dataset

   colnames(penguins)

   [1] "species"           "island"            "bill_length_mm"   
   [4] "bill_depth_mm"     "flipper_length_mm" "body_mass_g"      
   [7] "sex"               "year"             

2. To access any column variable, we use the $ syntax

   penguins$species[1:10]

    [1] Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie Adelie
   Levels: Adelie Chinstrap Gentoo

3. To get a table count of a variable or a cross-table count of 2 variables

   table(penguins$sex, useNA = "ifany")

   
   female   male   <NA> 
      165    168     11 

4. See Data Cleaning for more ways to work with dataframes

Data structures (dataframe, loading data)

1. To import data from csv file, use read_csv
2. To import data from an excel file, use read_excel
   
3. To import data from SAS/SPSS/Stata or export from R to these file formats, use the haven package

Data structures (dataframe, adding labels)

1. To add labels to a dataframe, we make use of the labelled package

Functions

1. A good chapter on explaining functions is Chapter 6, Advanced R, where most of the examples in the slides are taken

2. Example function

   f01 <- function(x, y) {
     # A comment
     x + y
   }

3. Components of a function

   1. Arguments \(\rightarrow\) inputs to the function (there can also be functions with no arguments)
   2. Body \(\rightarrow\) the code inside the function
   3. Environment \(\rightarrow\) a list-like structure that stores/maps the variable name to value (a namespace as some people call it)

4. Functions are objects that can be assigned or they can be anonymous functions

Function invocation

1. How do we invoke the function?

   print(f01)

   function(x, y) {
     # A comment
     x + y
   }

   f01(3, 5)

   [1] 8

2. How do we invoke one function after another?

   1. Suppose we have another function f02

      f02 <- function(a, b) {
        a * 10 + b
      }

   2. To call f01 first followed by f02 and sqrt we will “nest” the functions (inside-out, right to left)

      sqrt(f02(f01(1, 2), 5))

      [1] 5.91608

Function invocation (continued)

Another way to invoke a series of functions is to use the pipe operator %>% provided by the magrittr package

library(magrittr)
value <- 1
value %>%
  f01(2) %>%
  f02(5) %>%
  sqrt()

[1] 5.91608

Functions - argument passing

1. Pass by copy vs pass by reference

2. In R, all functions are pass by reference if you don’t modify the object subsequently within the function. i.e. copy-on-modify

3. How to pass arguments?

   f04 <- function(x, y, z) {
     (x / y) * z
   }

   1. By position e.g. f04(10, 50, 20)

   2. By name e.g. f04(z = 20, y = 50, x = 10)

   3. Unpacking multiple named arguments using do.call

      do.call(f04, list(x = 10, y = 50, z = 20))

      [1] 4

Common functions - dplyr::select

The select function enables you to select a subset of the columns in your dataset

# library(tidyverse) # if not already loaded
# library(penguins) # if not already loaded
penguins %>%
  select(species, island, bill_length_mm) %>%
  head(3)

# A tibble: 3 × 3
  species island    bill_length_mm
  <fct>   <fct>              <dbl>
1 Adelie  Torgersen           39.1
2 Adelie  Torgersen           39.5
3 Adelie  Torgersen           40.3

# there is also ends_with
penguins %>%
  select(starts_with("bill")) %>%
  head(3)

# A tibble: 3 × 2
  bill_length_mm bill_depth_mm
           <dbl>         <dbl>
1           39.1          18.7
2           39.5          17.4
3           40.3          18  

penguins %>%
  select(matches("mm")) %>%
  head(3)

# A tibble: 3 × 3
  bill_length_mm bill_depth_mm flipper_length_mm
           <dbl>         <dbl>             <int>
1           39.1          18.7               181
2           39.5          17.4               186
3           40.3          18                 195

Common functions - dplyr::filter

The filter function enables you to select a subset of the rows that meet a certain criteria

penguins %>%
  filter(body_mass_g > 3500) %>%
  head(3)

# A tibble: 3 × 8
  species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
  <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
1 Adelie  Torgersen           39.1          18.7               181        3750
2 Adelie  Torgersen           39.5          17.4               186        3800
3 Adelie  Torgersen           39.3          20.6               190        3650
# ℹ 2 more variables: sex <fct>, year <int>

penguins %>% filter(sex == "female")

# A tibble: 165 × 8
   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
 1 Adelie  Torgersen           39.5          17.4               186        3800
 2 Adelie  Torgersen           40.3          18                 195        3250
 3 Adelie  Torgersen           36.7          19.3               193        3450
 4 Adelie  Torgersen           38.9          17.8               181        3625
 5 Adelie  Torgersen           41.1          17.6               182        3200
 6 Adelie  Torgersen           36.6          17.8               185        3700
 7 Adelie  Torgersen           38.7          19                 195        3450
 8 Adelie  Torgersen           34.4          18.4               184        3325
 9 Adelie  Biscoe              37.8          18.3               174        3400
10 Adelie  Biscoe              35.9          19.2               189        3800
# ℹ 155 more rows
# ℹ 2 more variables: sex <fct>, year <int>

Common functions - dplyr::mutate

The mutate function enables you to add columns to your dataset. The added columns can be derived from existing column(s)

penguins %>%
  mutate(body_mass_100g = body_mass_g / 100) %>%
  select(body_mass_g, body_mass_100g) %>%
  head(5)

# A tibble: 5 × 2
  body_mass_g body_mass_100g
        <int>          <dbl>
1        3750           37.5
2        3800           38  
3        3250           32.5
4          NA           NA  
5        3450           34.5

penguins %>%
  mutate(bill_length_plus_depth_mm = bill_length_mm + bill_depth_mm) %>%
  select(matches("bill")) %>%
  head(5)

# A tibble: 5 × 3
  bill_length_mm bill_depth_mm bill_length_plus_depth_mm
           <dbl>         <dbl>                     <dbl>
1           39.1          18.7                      57.8
2           39.5          17.4                      56.9
3           40.3          18                        58.3
4           NA            NA                        NA  
5           36.7          19.3                      56  

Common functions - dplyr::group_by, summarise

Common functions - dplyr::if_else, case_when

The functions if_else and case_when are often used with mutate to create new variables

set.seed(2)
penguins %>%
  mutate(bill_length_type = if_else(bill_length_mm >= 48.5, "long bill", "short bill")) %>%
  sample_n(4)

# A tibble: 4 × 9
  species   island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
  <fct>     <fct>           <dbl>         <dbl>             <int>       <int>
1 Chinstrap Dream            43.5          18.1               202        3400
2 Gentoo    Biscoe           43.6          13.9               217        4900
3 Gentoo    Biscoe           48.1          15.1               209        5500
4 Gentoo    Biscoe           46.8          14.3               215        4850
# ℹ 3 more variables: sex <fct>, year <int>, bill_length_type <chr>

penguins %>%
  mutate(bill_length_type = case_when(bill_length_mm >= 48.5 ~ "long bill",
                                      bill_length_mm < 48.5 ~ "short bill",
                                      TRUE ~ NA_character_)) %>%
  select(bill_length_type) %>%
  table(useNA = "ifany")

bill_length_type
 long bill short bill       <NA> 
        87        255          2 

How are categorical variables represented?

1. Categorical variables are characters, but we need a numeric matrix in order to perform computations to get the \(\beta\) coefficients

set.seed(10)
idx <- sample(1:344, 5)
penguins[idx,]

# A tibble: 5 × 8
  species   island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
  <fct>     <fct>              <dbl>         <dbl>             <int>       <int>
1 Adelie    Dream               35.6          17.5               191        3175
2 Chinstrap Dream               50.7          19.7               203        4050
3 Adelie    Torgersen           39.7          18.4               190        3900
4 Gentoo    Biscoe              43.2          14.5               208        4450
5 Gentoo    Biscoe              47.2          13.7               214        4925
# ℹ 2 more variables: sex <fct>, year <int>

1. So one-hot encoding / dummy encoding will be applied to categorical variables

model.matrix(~species, penguins)[idx,]

    (Intercept) speciesChinstrap speciesGentoo
137           1                0             0
330           1                1             0
72            1                0             0
211           1                0             1
271           1                0             1

Framework

flowchart LR
  A[Model specification] --> B[Inference] --> C[Diagnostics]

1. Model specification
   1. How to decide what type of model to fit?
   2. How to decide what variables to choose?
2. Inference
   1. How do we get the value of the coefficients, confidence intervals, p-values?
3. Diagnostics
   1. Does the model fit our data well?

Why the need for GLM?

1. We want to find the association between an outcome (e.g. BMI) and some dependent variables (e.g. age, gender, ethnicity, social economic status)
2. We are used to fitting a multivariable linear regression \[ \begin{align} \underbrace{BMI}_{\text{continuous}} &= \underbrace{\beta_0 + \beta_1age + \beta_2gender + \beta_3ethnicity + \cdots + \epsilon}_{\text{this linear combination is continuous}} \\ \epsilon &\sim N(0, \sigma^2) \end{align} \]
3. What if your outcome is not continuous?
   1. E.g. outcome takes on values of 0 or 1 (logistic regression)
   2. E.g. outocme takes on discrete integer values 0, 1, 2, … (poisson regression)
4. GLM gives you a way to relate the outome (which can take on various distributions) with a linear combination of dependent variables

Family member 1: Linear regression

When our outcome, \(y_i\) is continuous and takes on real values (i.e. \(y_i \in \mathbb{R}\)), we may choose to model \(y_i\) with a normal distribution

Running a linear regression in R

1. Notice the use of the formula syntax. LHS of the ~ is the dependent variable, while RHS are the independent variables
2. Notice the family argument contains the 2 components of the GLM we discussed earlier

model <- 
  glm(body_mass_g ~ species + bill_length_mm +
                    bill_depth_mm + 
                    flipper_length_mm + sex,
      data = penguins,
      family = gaussian(link = "identity"))

summary(model)

Call:
glm(formula = body_mass_g ~ species + bill_length_mm + bill_depth_mm + 
    flipper_length_mm + sex, family = gaussian(link = "identity"), 
    data = penguins)

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -1460.995    571.308  -2.557 0.011002 *  
speciesChinstrap   -251.477     81.079  -3.102 0.002093 ** 
speciesGentoo      1014.627    129.561   7.831 6.85e-14 ***
bill_length_mm       18.204      7.106   2.562 0.010864 *  
bill_depth_mm        67.218     19.742   3.405 0.000745 ***
flipper_length_mm    15.950      2.910   5.482 8.44e-08 ***
sexmale             389.892     47.848   8.148 7.97e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 82563.33)

    Null deviance: 215259666  on 332  degrees of freedom
Residual deviance:  26915647  on 326  degrees of freedom
  (11 observations deleted due to missingness)
AIC: 4723.9

Number of Fisher Scoring iterations: 2

Family member 4: Conway-Maxell-Poisson (CMP) regression

1. When our outcome variable is a “count”, the poisson regression model is used because the poisson distribution models counts better than say linear regression which assumes that the outcome variable is normally distributed
2. However, the poisson regression model is restrictive as it assumes equi-dispersion in our outcome variable (i.e. mean of outcome variable approximately equals the variance of the outcome variable) When there is under- or over-dispersion, the poisson distribution no longer models the outcome variable well.
3. When there is over-dispersion (variance of outcome > mean of outcome), negative binomial regression is often used
4. When there is under-dispersion (i.e. variance of outcome < mean of outcome), we will use the Conway-Maxwell-Poisson regression model
5. Useful resources: (i) Video explanations see^1,2 (ii) Sellers & Premeaux³ explains the CMP regression and surveys available statistical packages

CMP distribution

1. The Conway-Maxwell-Poisson distribution differs from the Poisson distribution with the introduction of a parameter \(\nu\)
2. The \(\nu\) parameter enables the Conway-Maxwell-Poisson distribution to model a wide range of dispersion, and generalises well-known distributions
   1. When \(\nu = 1\), it takes on the Poisson distribution (equi-dispersion)
   2. When \(\nu = 0\) and \(\lambda < 1\), it takes on the geometric distribution (over-dispersion)
   3. When \(\nu \rightarrow \infty\), it takes on the bernoulli distribution (under-dispersion)

CMP regression

1. Use of the Conway-Maxwell-Poisson (CMP) distribution in the generalised linear model (GLM) setting:
   1. In the original CMP formulation of Sellers and Shmueli⁴, the relationship between dependent variable \(Y\) and independent variables \(X\) is given by the equations:
      1. \(ln(\lambda_{i}) = X_{i}^T\beta\)
      2. \(\mathbb{E}(Y_{i}) = \lambda_{i} \frac{\partial ln Z(lambda_{i}, \nu)}{\partial \lambda_i} \approx \lambda_{i}^{\frac{1}{\nu}} - \frac{\nu - 1}{2\nu}\), where \(Z(\lambda_i, \nu_i)\) is a normalizing constant
2. However, the CMP regression model is not amenable to interpretation of coefficients by mean contrasts, because \(\lambda_i\) is related to \(\mathbb{E}(Y_i)\) by a non-linear function
3. To overcome this, Huang 2017⁵ reparameterised the CMP distribution to model the mean directly. He called it CMP_\(\mu\)

CMP_\(\mu\)-regression

1. Under the CMP_\(\mu\)-regression formulation, the relationship between dependent variable \(Y\) and independent variables \(X\) is given by \(ln(\mathbb{E}(Y_i|X)) = X_{i}^T\beta\)
2. The convenient thing about CMP-regression and CMP_\(\mu\)-regression is that we can conduct a likelihood ratio test^4,5 to see if a poisson regression model (poisson regression is a special case of CMP-regression when \(\nu = 1\)) is adequate (i.e. \(H_0: \nu = 1\) vs \(H_1: \nu \ne 1\))
3. Note that \(H_1\) does not specify the direction (under vs over) of data dispersion. This can be assessed via the dispersion estimate \(\hat{\nu}\)⁴ (i.e. \(\hat{\nu} > 1\) if under-dispersion and \(\hat{\nu} < 1\) if over-dispersion)
4. Model diagnostics⁵: If the underlying distributional model is correct, the probability inverse transformation (PIT) should resemble a random sample from a standard uniform distribution. Goodness-of-fit can then be assessed by plotting a histogram of the PIT, or via a quantile plot of the PIT against the uniform distribution

Motivation: why do we care about missing data?

1. Commonly used statistical models (e.g. GLM) require complete-cases \(\rightarrow\) forced to discard data points
2. Data is expensive to collect
3. In throwing out data, you will face the problem of:
   1. Efficiency: number of data points are reduced \(\rightarrow\) lose statistical power
   2. Bias: you might systematically exclude a population group and end up with a biased dataset
   3. Both bias and efficiency will influence your final conclusion

Can these missingness mechanisms be tested for?

1. MAR and MNAR cannot be tested for (see section 1.4 of Enders¹⁶)
   1. In each case, an independence assertion regarding the missing variable was made
   2. Since the missing variable is by definition partially unobserved, we cannot do a test for such an assertion where the input to the test is unobserved
2. MCAR can be tested for (see section 1.9 of Enders¹⁷)
   1. There are testable propositions, the idea is to assume that the cases with missing data came from the same population and thus have the same means and covariances as the complete cases
   2. Little’s test is one test for MCAR that evaluates the mean differences across subgroups of cases that share the same missing data pattern
   3. The null hypothesis for Little’s test states that the data are MCAR, so a statistically significant test statistic provides evidence against the MCAR mechanism

Caveats about Little’s test

From section 1.9 of Enders¹⁷

1. The test does not identify the specific variables that violate MCAR, so it is only useful for testing an omnibus hypothesis that is unlikely to hold in the first place
2. The version of the test outlined above assumes that the missing data patterns share a common covariance matrix. MAR and MNAR mechanisms can produce missing data patterns with different variances and covariances, and the test statistic in Equation 1.4 would not necessarily detect covariance based deviations from MCAR
3. Simulation studies suggest that Little’s test suffers from low power, particularly when the number of variables that violate MCAR is small, the relationship between the data and missingness is weak, or the data are MNAR. Consequently, the test has a propensity to produce Type II errors and can lead to a false sense of security about the missing data mechanism

R Markdown

1. R Markdown is an extension to markdown that enables user to execute R code and display the code output in addition to text
2. The knitr application enables us to execute R code and display output in a temporary markdown document
3. Pandoc then converts the markdown document to the desired format (e.g. html, pdf, docx, pptx)
4. See summary by Robin Linacre

Image from R Markdown cookbook chapter 2.1

R markdown YAML header

---
title: "My project"
author: "CRU"
date: '`r format(Sys.Date(), "%d %b %Y")`'
output: 
  bookdown::html_document2:
    toc: true
    toc_float: true
bibliography: references.bib
csl: vancouver-superscript.csl 
link-citations: true
---

1. This is a common setting which I use
2. bookdown::html_document2 of the Bookdown package is a good output format because it takes care of formatting e.g. auto-numbering of sections, handling of references, appendix
3. toc (line 7 & 8) refers to table of contents
4. bibliography, csl and link-citations is for adding citations stored in a bibtex file into your document. More will be covered in a later slide